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We study non-axisymmetric oscillations of rapidly and differentially rotating relativistic stars in 
the Cowling approximation. Our equilibrium models are sequences of relativistic polytropes, where 
the differential rotation is described by the relativistic ji-constant law. We show that a small degree 
of differential rotation raises the critical rotation value for which the quadrupolar f-mode becomes 
prone to the CFS instability, while the critical value of T/|W| at the mass-shedding limit is raised 
even more. For stiffer equations of state these effects are even more pronounced. When increasing 
differential rotation further to a high degree, the neutral point of the CFS instability first reaches 
a local maximum and is lowered afterwards. For stars with a rather high compactness we find 
that for a large degree of differential rotation the absolute value of the critical T/| W] is below the 
corresponding value for rigid rotation. We conclude that the onset of the CFS instability is eased 
for a small degree of differential rotation and for a large degree at least in stars with a higher 
compactness. Moreover, we were able to extract the eigenfrequencies and the eigenfunctions of r- 
modes for differentially rotating stars and our simulations show a good qualitative agreement with 
previous Newtonian results. 

PACS numbers: 04.30.Db, 04.40.Dg, 95.30.Sf, 97.10. Sj 

I. INTRODUCTION 

Rapidly spinning neutron stars can be destabilized due to the CFS [TH3] instability. This instability, for quadrupole 
non-axisymmetric perturbations with m = 2, sets in for uniformly rotating polytropes when (3 = T/|W| ~ 0.14 where 
T is the rotational kinetic energy and W is the gravitational potential energy of the star. For higher values of m the 
instability sets in for smaller rotational rates but the growth time is considerable longer and viscosity stabilizes the 
star. The CFS instability is generic for r-modes, i. e. these modes are CFS unstable for any rotational rate while for 
the f- modes the previous value of /3 is approximately correct. This value of j3 corresponds to rotational frequencies 
of the stars which are 85-90% of the Kepler frequency. These high spin frequencies have not yet been observed in 
nature but it is expected that they may be met in newly born neutron stars. Newly born neutron stars are expected 
to rotate not only fast but also differentially, and this is a factor that has to be taken into account in the study of 
instabilities. 

Actually, differential rotation can appear in many phases of the stellar evolution of a neutron star, such as in 
protoneutron stars |H [S], in the massive remnant of binary neutron star mergers [SHE], or as a results of stellar 
oscillations (r-modes) that may drive the star into differential rotation via nonliner effects [9TTT2"]. The phase where 
the star rotates differentially can last between seconds to months, depending on the dissipative mechanism that drive 
the star to uniform rotation, such as viscosity, magnetic braking |T31 H3] and turbulent motion [T5J [TC]. However, 
this very short period is met during the most violent phases of neutron star's life, as in the case of a core collapse 
or binary merging. This is exactly the time when we expect to get the strongest emission in gravitational waves. 
Since the ground based detectors are reaching sensitivities which allow the detection of gravitational wave signals from 
oscillating or unstable neutron stars the critical point for the onset of instabilities, the growth times of the instabilities 
and the exact frequencies of the emitted waves are urgently needed. 

In this work we study within the general relativistic framework the oscillations and instabilities of fast and differ- 
entially rotating neutron stars in the so called Cowling approximation i.e. assuming a fixed spacetime. This is the 
first study of its kind since earlier works were using certain approximations e.g Newtonian theory [T7] or slow rotation 
[TUE3. Results, from fully nonlinear calculations exist only for the axisymmetric oscillation [201 HI]- Finally, the 
f-mode and the secular stability limits have been investigated in |22j . 

Recently, it has been found that in stars that rotate with a high degree of differential rotation an m — 2 dynamical 
instability can appear even for considerably low rotation rates (T/\ W\ ~ O(10 -2 )) as suggested in [231I21]- In addition, 
an m — 1 dynamical instability has been identified for high degrees of differential rotation and soft equations of state 
[25l [26] . More recently, the discovery of m — 1 and m — 2 dynamical instability even for stiff equations of state [27] 
has been reported . Studies based on linear analysis [35] suggest that low T/|W| instabilities might be triggered 
when the corotation points of the unstable modes fall within the differentially rotating structure of the star. 

The structure of the paper is as follows. In Section [Tl] we present the formulation that we have used to construct 
the equilibrium configurations. We then derive the perturbation equations in the general case of barotropic and 
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non-axisymmetric perturbations, which are numerically solved in Section III for the non-axisymmetric and barotropic 
case. Finally, in Section IV we summarize the crucial results of this study. 
Throughout the paper we use geometrical units c = G = 1. 



II. FORMULATION OF THE PROBLEM 



In this section we first show the basic principles of our problem. We then derive the evolution equations governing 
oscillations of perturbations on compact objects and the appropriate boundary conditions. Next, we present the 
equations of state we use for constructing equilibrium models and what kind of sequences we will consider. Finally, 
we explain how we solve the problem numerically. 



A. Basic Principles 



We study linear perturbations of rapidly and differentially rotating, relativistic neutron stars on a fixed background. 
First of all, this requires the computation of equilibrium configurations which are prescribed by their general-relativistic 
line-element and the specific form of the energy-momentum tensor. This last quantity describes the proper charac- 
teristics of the neutron star fluid and therefore, a certain equation of state has to be specified as well; we will address 
this issue later in this section. 

The background metric of a fast rotating, relativistic star has to be stationary and axisymmetric and using spherical 
coordinates (r,0,<p), the line-element can be written as 

ds 2 = -e 2v dt 2 + e 2 V sin 2 6(d(f> - udt) 2 + e 2 ^{dr 2 + r 2 d9 2 ) . (1) 

The four unknown functions ip, fx, v and uj are the metric potentials which depend on r and 6 only and can be 
determined by solving the time-independent Einstein equations. Furthermore, for equilibrium models the radial and 
the polar component of the fluid's four-velocity vanishes while the remaining two components are related by the 
angular velocity fi via 

(u\ u r , u e , u*) = (u\ 0, 0, flu*) . (2) 

In general, the angular velocity is a function of r and 6 and although the exact form of r2(r, 9) is arbitrary to 
some extent, we adopt the so-called relativistic j-constant rotation law 

[ c ' ~ 1 - (fi - Lofr^ sin 2 fleW"") ' ( ' 

which satisfies the Rayleigh criterion for local dynamical stability against axisymmetric perturbations and which is 
commonly used in similar studies. [5DJ |3D1 131] • Here, f2 c is the angular velocity on the rotation axis and the free 
parameter A controls the degree of differential rotation. This quantity has units of length in the Si-system and 
specifies the length scale over which the angular velocity varies inside the star. As one can see from equation ([3]), the 
star is rotating uniformly again when A tends to infinity. 

This definition of A is not convenient when considering sequences of equilibrium models with increasing rotation 
rate, since the radius of the star can vary by a factor of two which would be reflected in different degrees of differential 
rotation for each model. This issue can be circumvented by using a slightly modified parameter 

A := A/r e , (4) 

which normalizes A by the equatorial radius r e of the star; see also |20j . This definition maintains the same degree of 
differential rotation for all models along a sequence. 

We furthermore assume the neutron star to be a perfect fluid without viscosity. In this case, the energy-momentum 
tensor has the form 

= (e + p)vfu v +psT , (5) 

where e is the energy density, p is the pressure and is the inverse metric. Energy and pressure are not independent 
quantities but they are related by an equation of state. In this study, we will focus on polytropic equations of state 
in the form of 

P = Kp T , (6) 
e = p + Np . 



3 



Here, p is the rest-mass density, K is the polytropic constant, N is the polytropic index and T = 1 + 1/N. For 
barotropic oscillations, pressure- and density-perturbations are connected by the speed of sound; Sp = c 2 Se. In the 
case of polytropic equations of state, its value can be computed analytically and it is 

ol = (7) 
s e + p v ' 

The thorough study of neutron star oscillations is very demanding and requires to solve the full, non-linear equations 
of general relativity. At this point, we will introduce some approximations in order to deal with a smaller set of 
equations and unknowns. First of all, we restrict the study to small perturbations around the equilibrium which 
allows us to linearize the equations. Second, we will use the so-called relativistic Cowling approximation, in which all 
metric perturbations are neglected. This effectively discards the space-time degrees of freedom and one only has to 
deal with the equations governing the motion of the fluid, which follow from the linear variation of the local law of 
energy- momentum- conservat ion 

8(V V T^)=0. (8) 

Since in the Cowling approximation, the covariant derivative is not affected by the fluid perturbation, equation ^ 
can be written as 

dJT^ + T^JT 1 " + T v yv ST^ = , (9) 
where ST 11 " is the perturbed energy-momentum tensor. 



B. Perturbations Equations 



The perturbation equations can be derived from equation (|9) and for this, we will follow a very convenient formu- 
lation developed in 32, 33J. Instead of using the canonical fluid perturbation variables like velocity- or pressure- 
perturbation, the independent components of the perturbed energy-momentum tensor are utilized for the time- 
evolution. 

To be more specific, we are directly evolving the variables Q\ := ST lt , Q 2 := ST*^, Q3 := ST tr and Q4 :— ST in 
time, using equation ([9| to obtain the corresponding evolution equations. Additionally, we make use of an auxiliary 
variable Qe ■= 5T rr which turns out to be a linear combination of Q\ and Q§ and which helps us in writing the 
resulting set of differential equations in a more compact way. 

The final set of evolution equations is then given by 



d t Qi = - d rj ,Q 2 - d r Q 3 - deQi 



= -2i/+2t/> 2 



r sin 9(il — uj)d r Lu + d r ip + 2<9 r /i + ?>d r v + 



= _ 2l ,+2i/> 2 



r 2 sin 2 6(0, - uj)d e uj + d e i> + 2d e n + 3d e v 



/•_ 
cos f 
sin 9 



(10a) 



d t Q 2 = 2 d 4 ,Q 1 ~ 208^2 - Od r Q 3 - OdeQi + 



e- 2 »+ 2 »(0-u>) 2 - 



— 2v+2ib 2 • 2 



+20d e p + (0 + 2uj)d e v + (SO - 2uj)- 



cos 9 



sin 9 



r 2 sin 2 9 



r 2 sin 2 9(0 - u)ujd r uj + d r (O - u) + (30 - 2oj)d r ip 



-20d r p +(0 + 2uj)d r v + (20 - w) 



D -2»+24> r 2 sin 2 e( p _ ^deu + d g (0-uj) + (30 - 2uj)d 9 i> 



(10b) 
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dtQz = ~ Sld^Qz - d r Q 6 

- [e- 2 ^ +2u d r u + e^ +2 ^rsin 2 9 [(SI 2 - uj 2 ) (1 + rd r ij>) - ru>d r ui]] Q t 
_ [ e -2M+2V V sin 2 e [ rdrU _ 2 (n _ + rd r ip)]] Q 2 

d r u + 2d r n + e^+^r sin 2 (SI - u) 2 (1 + rd r ipj\ Q 6 



(10c) 



Sld^Qi - ^d e Q 6 



[ e -' 2 ^ +2v 'd 6 v + e - 2 M+2V> r 2 sin6) |-^2 _ w 2^ (cos^ + sin^eV 7 ) -smeudeu]] 
[ e -^+^r 2 sin 6» [sin 9d oj - 2 (ft - w) {cos0 + sm6d e ip)]] Q 2 
d B v + 2d 6 [L + e ~ 2v+2 ^r 2 sin 9 (Q - uj) 2 (cos 9 + sin 6d g ip) 



(lOd) 



Additionally, the value of Q e in terms of Qi and Q 2 is calculated by 



1<A = 



1 ~ e -2^+2^ r 2 sin 2 ,9(Q _ w ) 



2„2 



(11) 



1 



-2i/+2V„2 



sin 2 6»(fi 2 - uj 2 )} Qi - 2eT 2vJr2 ^r 2 sin 2 9 (SI - oj)Q 2 ] 



Since equilibrium stars are axisymmetric, but not neccessarily spherically symmetric in case of rapid rotation, we 
will not use an angular decomposition into spherical harmonics but instead only separate the azimuthal part according 
to 

Q i (t,r,6,<t>):=Q i (t,r,9)e im ' t > for i=l,...,4. 



Due to this decomposition, the evolution equations ( 10a |-( lOd) are only slightly modified and we just have to perform 



the substitution — > im. In the following discussion, we will omit the tilde to avoid confusion. 



C. Boundary Conditions 



In order to close the system of equations (10a)-(10d), one needs to impose proper boundary conditions. In the 



study presented here, this boundary consists of the stellar surface, the equatorial plane and the rotation axis. Special 
care has to be taken of the origin which is mapped onto a boundary line in spherical coordinates. 

Let us first consider the stellar surface, where the radial and polar perturbation variables Q3 and Q4 are zero by 
virtue of their definition 

Q 3 = ST tr = (e+p)u t Su r = 
Q 4 = ST te = (e + p) u'Su 9 = , 

because pressure and energy density vanish there. Qi and Q 2 however have to remain continuously differentiable 
there. 

On an axisymmetric background, each oscillation mode can be assigned to one of two possible parity classes which 
behave differently under space inversion; see [34]. Due to the decomposition with respect to the azimuthal angle <f>, 
the behaviour under reflection with respect to the equatorial plane is also well-defined and the two parity classes obey 
different conditions in the equatorial plane. For one class of modes, the perturbation variables Q± 7 Q 2 and Q3 are all 
even, whereas Q4 is odd. The I = m fundamental mode belongs to this class, since the angular part of the pressure 
perturbation behaves like the scalar spherical harmonic Yj m , which is even as well. In contrast, the perturbation 
variables of modes belonging to the second class behave exactly the other way around. The / = m r-mode is an 
example for a mode of this specific class. We therefore apply different boundary conditions along the equatorial plane, 
depending on which class of modes we want to excite. 

For deriving the boundary conditions along the rotation axis and at the origin, we use a representation of pressure 
and velocity perturbations according to |35j . These expressions were derived in the so-called slow-rotation approxi- 
mation where one neglects high order effects of rotation, but the relevant coupling between the two different mode 
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classes is already included within this approach and does not change at all in rapidly rotating stars. The bound- 
ary conditions we use in this study for discussing spherically symmetric (mostly for code verification purposes) and 
quadrupolar perturbations are shown in Tables |T| and |TT) 



finite finite 

2 



TABLE I: Boundary conditions for the perturbation variables at the origin. 



finite finite finite 

2 d e =0 d e = 



TABLE II: Boundary conditions for the perturbation variables along the rotation axis. 



D. Equilibrium Models 

For constructing the equilibriums models, we use the rns-code developed by Stergioulas |36l [37] , which was later 
extended to handle differential rotation. Furthermore, we are using three different equations of state which were 
already used for the study of uniformly rotating stars in [38, 39 and for which we construct equilibrium sequences 
with a constant central rest-mass density. All three EoS are polytropic ones; but two of them are fits to tabulated 
equations of state and together they cover very well the expected range of neutron star massses and radii, see for 
example [4TJIl4"T] . 

For all equations of state, we consider uniformly rotating models as well as their differentially rotating counterparts 
for different values of A. All sequences start with a non-rotating model for which the axis ratio r p /r e of polar to 
equatorial coordinate radius equals 1. Subsequent models are computed by decreasing this ratio by a factor of 0.05 
until the Kepler limit is reached. 

The EoS B is a widely used polytropic equation of state with parameters T = 2 and K = 100 in units of G = 
c = Mq = 1, and here we consider the commonly used BU sequence, for which the dimensionless central rest-mass 
density is given by p c = 1.28 x 10~ 3 . This leads to standard values for the mass and the radius of the non-rotating 
model; a detailed list with properties of the BU models and their differentially rotating couterparts with A = 1 (B 
sequence) is given in [20]. The non-rotating model, which is labelled BU0, has a gravitational mass of M = 1.4 Mq 
and a circumferential radius of R e = 14.16 km. The fastest rotating member of the uniformly rotating B sequence has 
an axis ratio of r p /r e = 0.58, a slightly increased mass of M = 1.695 Mq and an equatorial radius of R = 19.96 km. 

When the star is differentially rotating (the following discussion refers to the B sequence with A = 1), these 
changes become even more evident. Due to the fact that differential rotation allows to store a large amount of angular 
momentum near the rotation axis, the value of T/|VT| can reach much higher values. 

Furthermore, the star can become more compact near the rotation axis in comparison to the corresponding uniformly 
rotating star which results in a very high mass. The axis ratio for the most rapidly rotating model B12 can be decreased 
down to a value of 0.4 with a mass of M = 2.532 M which is an increase of 50 % when compared to the uniformly 
rotating model BU9. 

As already mentioned, we additionally study two other equations of state which are polytropic fits to the tabulated 
EoS A and EoS II, see 42444] for more details. The oscillation frequencies of neutron star models governed by these 
EoS have already been investigated in [35]. The fitting parameters for EoS II are T = 2.34 and K = 1186 and T = 2.46 
and K — 1528 for EoS A again in units where G = c = Mq = 1. For both EoS, we choose models which are very 
close to their maximum mass; e.g. see Figure 9 in [3S]. Apparently, the models for EoS A and EoS II have smaller 
radii and larger masses compared to the B sequences. A detailed list of the differentially rotating background models 
is shown in Table ITITl 

E. Numerical Implementation 

After decomposing the perturbed quantities with respect to cj), the resulting evolution equations form a two- 
dimensional problem in the spherical coordinates r and 9. As radial grid coordinate we choose a mapping of the 
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TABLE III: Stellar properties of the differentially rotating models with A = 1. f2 c and fi e are the angular velocities along the 
rotation axis and at the equator respectively, R e is the equatorial circumferential radius and M the gravitational mass. 

form s = 0.5 (r/r e ) 2 , which allows for stable time evolutions with much less artifical viscosity when compared to the 
original radial grid coordinate of the rns-code. For the angular grid coordinate we choose between t — cos(0) and 
t = 1 — 29 /n depending on the degree of differential rotation. In the Cowling approximation, only the interior of the 
star needs to be considered and this means that s < 0.5 for all t. In the non-rotating limit, the surface of the star 
coincides exactly with the gridline at s = 0.5. 

However, this choice of grid parameters poses a small challenge when applying boundary conditions on the surface 
of a rapidly rotating star, since in this case the surface lies in between the grid points of our computational domain. 
One therefore has to approximate the outer boundary in a proper way which is done here in the following manner: 
For each angular direction, we monitor the variation of the energy density when moving radially outwards and set the 
surface at this particular angle to be the first grid point where the energy density vanishes. The boundary conditions 
are then applied to these set of grid points. Depending on the shape of the star, one may have to increase the angular 
resolution of the grid in order to get an accurate approximation of the stellar surface. As already mentioned before, 
we use two different angular grid coordinates. While the mapping t — cos(#) is appropriate for small degrees of 
differential rotation, for higher degrees t = 1 — 2 9/tt is the better choice. 

Concerning the time-evolution scheme, we use standard central differences for spatial discretization and the classical 
fourth-order Runge-Kutta method as time integrator. Similar to previous studies within the same perturbative 
framework [38, 3ET, this code is prone to numerical instabilities around the origin of the star, which can be cured by 
using standard second-order Kreiss-Oliger dissipation 45 . Depending on the resolution of the grid, we always use 
the smallest amount of artificial viscosity in order to keep the effect of this purely numerical additive as marginal as 
possible. Typically, the simulations were performed using a grid size of 120 x 90 points; in this case, the dissipation 
coefficients are of the order of 10~ 4 . As it turned out, the mode frequencies do not change significantly when increasing 
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either the resolution of the computational domain or the amount of artificial viscosity. Additionally, the utilization of 
the Runge-Kutta method as time-integrator allowed for a further decrease in the amount of dissipation when compared 
to the Iterated Crank-Nicholson scheme used in previous studies. This also leads to a considerably slower numerical 
damping of the oscillations. 

The code is numerically stable for simulations up to some hundreds of milliseconds. Usually, our simulations cover 
roughly 100 ms, resulting in a frequency resolution of A/ s» 10 Hz. We then extract the eigenfrequencies from the 
time evolution by perfoming a Fast Fourier Transform on the obtained time series at some arbitrary points inside the 
star. For an unambiguous identification of a specific mode, we use Fourier transforms at each single point inside the 
star in order to construct the entire two-dimensional eigenfunction. Once we determined an eigenfrequency by means 
of the Fourier spectrum at a certain point, we track the amplitude of the Fourier transform for all grid points at this 
particular frequency which gives us the modulus of the corresponding eigenfunction. A sign flip in the eigenfuction is 
reflected by a shift of tt in the Fourier transform's argument. With the eigenfunctions extracted, we are also able to 
perform mode recycling 21, 38J in order to enhance a particular mode in the oscillation pattern, typically resulting 
in a more accurate frequency determination. 



III. RESULTS 



We now present the results obtained with our method as outlined in the previous sections. First, we will test 
it on uniformly rotating models to check the accuracy of the newly developed code. We then turn to differentially 
rotating stars and investigate the effects for both small and large degrees of differential rotation on the onset of the 
CFS unstable fundamental quadrupolar mode. Finally, we will take a look at the generically unstable r-mode in 
differentially rotating stars. 



A. Uniform Rotation 



We first want to consider perturbations of uniformly rotating stars mostly for code verification purposes. Barotropic 
oscillations of this type have already been studied extensively in |38j with the very same polytropic equations of state 
considered in this study. However, here we use a completely different set of perturbation equations, different coordinate 
systems and numerical methods. Still, both axisymmetric and non-axisymmetric mode frequencies can be reproduced 
with a very good accuracy; the discrepancy between the two different codes is below 2% in all cases. 

As an example, Figure [l] shows a comparison between results obtained in the present study and published values 
for the fundamental quadrupolar f-mode splitting found in |38j . The sequences of background models in these two 
studies are identical, so a good means in order to estimate the differences between the two methods is to observe 
the off-centering of the triangles, which are the data points from the current study, from its enclosing circles, which 
represent the literature values. 



g~o Gaertig, Kokkotas (2008) 
"Kruegeretal.(2009) 




spin frequency / kHz 



FIG. 1: Comparison of the fundamental mode splitting for I = \m\ = 2 and different EoS. Differences can be monitored by the 
shift of the triangles (present study) with respect to its enclosing circles (results from |38|). 
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Obviously, only for very high rotation rates there are some deviations; otherwise the two codes provide nearly 
identical results. 

We also checked the frequencies of the second class of oscillation modes which are present in a barotropic, perfect 
fluid star as soon as rotation sets in; these are the inertial modes. Again, the agreement between the two codes is 
excellent; a comparison with frequencies found in [3D] also shows a good agreement (see [17] for more details). 



B. Axisymmetric Modes of Differentially Rotating Stars 



Next, we consider axisymmetric perturbations on differentially rotating backgrounds and compare the results with 
20J, who investigated the same problem with a non-linear code applied to small perturbations. 

Figure [2] shows the effect of rotation on the three fundamental modes F, 2 f and 4 f as well as the first overtones Hi 
and 2 pi for both uniformly and differentially rotating models of the BU and B sequences. As already discussed in 
Section [IIP I differential rotation allows for a higher mass of equilibrium models compared to their uniformly rotating 
counterparts and the corresponding value of T/|W| can reach significantly higher values. This is the reason why the 
dashed curves for uniformly rotating stars terminate already at a rather low T] \W\ in Figure [2j which is similar to 
Figure 5 in [3D]. 




0.00 0.05 0.10 0.15 0.20 

T/IWI 



FIG. 2: Mode frequencies of several quasi-radial (m = 0) modes for uniformly and differentially rotating stars of the B sequence. 



As one can clearly see, the mode frequencies of differentially rotating models are typically larger than the cor- 
responding frequencies for the uniformly rotating sequence. The reason for this is that differentially rotating stars 
usually have a smaller equatorial radius than their uniformly rotating counterparts with the same T j \ W\. Therefore, 
the crossing time of acoustic waves, like the f- and higher order p-modes, is lowered which leads to a higher oscillation 
frequency. 

The agreement with the data provided in [2D] is very good; the relative difference is below 2.5% for all frequencies. 
However, we disagree about the interpretation of the so-called Fn-mode which was found in [20] and was interpreted 
as a possible splitting of the fundamental quasi-radial mode in differentially rotating stars; a reason for this might be 
the Cowling approximation which violates energy and momentum constraints. 

When we reiterated this simulation with the code presented in this study, we also found peaks in the Fourier 
spectrum at the very same frequencies published in [2D] for the Fn-mode. However, after extracting the corresponding 
eigenfunction, we are convinced that this Fn-mode is nothing more than the 4 f-mode. The frequencies of both the 
F-mode and the f-mode are nearly indistinguishable in slowly rotating stars; the difference is around 100 Hz which 
makes a reliable assignment of the eigenfunctions a difficult task. Fortunately, it turned out that the eigenfunction of 
5u is a very good means to discriminate between the two quasi-radial modes. 

Figure [3] shows contour lines of the time-evolution variable Q4 ~ Su e (see section II C) projected on the (s,i)-plane 



for the F- and the f-mode. In both panels, these modes are compared for a very slowly rotating model (solid line) 
with an axis ratio of r p /r e = 0.999 and a B3 neutron star (dashed line). The step-like curve on the right side of both 
panels is the approximated surface of the B3 model. As already discussed in Section II E only for non-rotating stars 
the stellar surface coincides with s — 0.5. Within the grid-resolution used in Figure |3j this is still true for the slowly 
rotating model considered here. The rotation rate of the B3 model is also the reason, why the contours of Q4 are 
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shifted to smaller radial distances in this case. Apart from that, the shapes of the eigenfunctions are identical in both 
panels. 



rotation axis (t=l ) rotation axis (t=l) 




equatorial plane (t=0) equatorial plane (t = 0) 



FIG. 3: Left Panel: Contours of the Q4-eigenfunction of the F-mode for a slowly rotating B model with r p /r e = 0.999 (solid) 
and a B3 neutron star (dashed). The step-like curve on the right approximates the stellar surface for the B3 model. Right 
Panel: Same for the 4 f-mode. 

Actually, this issue about erroneous mode identification will always be present for various stellar models and certain 
modes. In the case shown here, the problem can be resolved if one drops the Cowling-approximation because this 
approximation overestimates the frequencies of the F-modes, making them nearly identical to those of the 4 /-mode 
for this particular sequence of equilibrium models. However, the problem will still be there for different modes or 
stellar models and thus the analysis presented here is the only meaningful way to distinguish between certain modes. 



C. Non-Axisymmetric f-modes in Differentially Rotating Stars 

We next turn to non-axisymmetric perturbations on differentially rotating stars. In the following discussion, we will 
mostly focus on the quadrupolar fundamental mode in particular, because this is the lowest order mode which is able 
to emit gravitational waves. For uniformly rotating stars, the 2 f-mode is also prone to the secular CFS instability, 
but typically only for models which rotate close to their mass-shedding limit. The neutral point for the onset of the 
instability is considerably lower for modes with higher azimuthal index m. However, it turns out that the growth 
time of these modes is considerably larger than the viscous timescales. Therefore, one usually studies the 2 f-mode 
exclusively. 



In Table IV we show the eigenfrequencies of the f-mode for the sequence B and a degree of differential rotation of 



A = 1.0. As one can see, the star becomes unstable at a point between the B8 and the B9 model. The values in Table 



IV have been compared with preliminary results from non-linear evolution codes and show a good agreement [48] . 

One may now ask, how this behaviour changes for different degrees of differential rotation. As we know, the 
uniformly rotating model, which corresponds to A~ x = 0.0, becomes secularly unstable just around its mass-shedding 
limit while for Ar 1 = 1.0 this happens at roughly 65% of the highest possible T/|W| for a stable equilibrium. 

Figure |4] shows the splitting of 2 f-mode for four different degrees of differential rotation, starting from the uniform 
rotation limit with A^ 1 = 0.0 to a value of A" 1 = 1.43. First of all, one can see that differential rotation leads 
to an increase in the mode frequencies compared to the uniformly rotating case; something we already observed 
for axisymmetric modes, see Figure [2j Second, a high degree of differential rotation favours the onset of the CFS 
instability in the sense that the neutral point is reached at a smaller fraction of the maximum allowed T/|W|. 

Third, for sufficiently large values of A -1 , the mode frequencies of the corotating branch always increase as can be 
seen from Figure [4] for A^ 1 > 1.0. This is in strong contrast to the uniformly rotating case where the frequencies of 
the prograde travelling modes reach a local maximum before they decrease again; see also Figure [2j If the background 
models allowed for high rotation rates, the frequency of the corotating branch would decrease even more and at some 
value of T/|W| merge with the frequency of the retrograde travelling mode which is dragged forward by rotation in 
the inertial frame. This point then signals the onset of the dynamical bar-mode instability. 
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TABLE IV: Frequencies of the 2 f-mode for the B sequence and A = 1.0. 



0.45 




FIG. 4: Left Panel: Splitting of the 2 f-mode for different degrees of differential rotation A. All equilibrium sequences are based 
on EoS B. Right Panel: Density contour lines of the most rapidly rotating model B12. The outermost line corresponds to the 
surface of the star, the z-axis represents the rotation axis. 

The situation changes in the case of differential rotation. Depending on the degree of differential rotation and the 
rotation rate, there is a transition of the axisymmetric equilibrium configuration from an oblate spheroid with its 
maximum energy-density and pressure at the center of the star towards a toroidal configuration where these peak 
values are reached along a ring-like structure in the equatorial plane; see the right panel of Figure [4] for a density- 
contour plot of the most rapidly rotating model with A^ 1 = 1.0. This change in the topology of the background 
model leads to an increase in the frequencies of the corotating branch which of course sets in earlier for equilibrium 
configurations with large degrees of differential rotation since in this case, the toroidal structure of the background star 
is already realized for comparatively small values of T/|W|. One should also note that while the sequences considered 
here are still configurations with the same central energy density, the maximum value of the energy density is reached 
outside the center of the star as soon as the transition to a toroidal structure sets in. 



1. Small Degree of Differential Rotation 

We want to discuss the onset of the CFS-instability a bit more in the regime of weakly differential rotation, i. c. 
small values of A -1 . As one can see clearly from Figure [ij the critical value of /3 := !T/|W|, where the neutral point 
of the secular instability is reached, typically increases with the degree of differential rotation. In order to estimate 
whether differential rotation favours the onset of the CFS mechanism, one also needs to know the maximum value /3 
can attain. In differentially rotating stars, a large amount of angular momentum can be stored in inner layers close 
to the rotational axis. This also means, that the highest possible value of j3 increases as well. 
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In the subsequent discussion, we will restrict our simulations to smaller values of A^ 1 for the following reason: 
Since we later want to normalize the critical value for the onset of the CFS instability j3 c with the corresponding value 
j3 s at the mass-shedding limit, we need dynamically stable equilibrium models at this Kepler frequency; otherwise 
the normalization would not make much sense. It turns out, that our particular choice of background models is 
dynamically unstable to radial oscillations for comparatively small values of A -1 . For example, the B sequence is 
dynamically unstable before reaching the mass-shedding limit for A^ 1 > 0.7. It is worth pointing out, that the 
endpoint of the B model series in Table |Hl} i.e. model B12, does not terminate because of the mass-shedding limit 
but because of dynamical instabilities. The same is also true for sequences A and II; due to the higher initial masses, 
this effect is even more pronounced there. 

Figure [5] shows, how /3 C and j3 s change with the degree of differential rotation. As already conjectured, both values 
increase with A~ l and the slope of j3 s is even larger than the one for j3 c . This confirms our initial assumption, that 
the onset of the secular CFS instability is eased for larger degrees of differential rotation. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 



A- 1 

FIG. 5: The values of /3 C and /3 S as a function of A -1 for the 2 f-mode and the B model series. The filled area marks the 
CFS-unstable region. 

In FigurejH] we plot the value of f3 c :— /3 C / f3 s for three different modes of the B sequence as function of the degree of 
differential rotation A -1 . This normalization of the critical T/|W| is similar to the common procedure for uniformly 
rotating stars, where one typically measures the angular velocity of the star in units of the Kepler-frequency. 



P 

1 c 




A" 1 

FIG. 6: The normalized, critical value /3 C for the onset of the CFS instability of the B sequences as a function of the degree of 
differential rotation A -1 for m = 2, 3, 4. 

Keep in mind that the uniformly rotating limit is reached for A^ 1 = 0.0. As we already have seen in Figure [4J in 
this case the 2 f-mode becomes unstable just around its mass-shedding limit which corresponds to j3 c f=a 1.0. It has 
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already been mentioned in Section III C that higher order modes become secularly unstable even earlier and this is 
also reflected in Figure [6] For example, in the limit of uniform rotation, the 3 /-mode gets unstable at roughly 70% 
of the critical T/\W\, the 4 f-modc even at approximately 50%. However, it has already been noted that the growth 
time for the instability will increase for modes with larger m and is most likely suppressed by viscosity. 

When we trace the value of /3 C for the various fundamental modes and several degrees of differential rotation, we 
see clearly that for larger values of A -1 , the critical T/|TT"| decreases invariably for all modes: i.e. the onset of the 
instability is favoured in these cases. 

We also studied the behaviour of (3 C for the stiffer equations of state EoS A and EoS II. The actual equilibrium 
models we considered are very close to their maximum allowed mass and therefore we can increase the differential 
rotation parameter A only to a very moderate degree of around A 2.3 because models with a higher degree of 
differential rotation will be dynamically unstable. Figure [7] again shows a plot of /3 C for different values of A -1 , 
but this time for the fundamental quadrupolar mode and the three different polytropic equations of state that were 
considered in this study. As already found in [35] and one can also see in Figure [l] the CFS instability acts much 
earlier in the more compact sequences A and II. 
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FIG. 7: The normalized, critical value /3 C for the onset of the CFS instability of the f-mode for different equations of state. 

It is worth noticing that the influence of differential rotation is even more pronounced for these sequences when 
compared to sequence B. As one can see from Figure [7J the critical curves for EoS A and EoS II decrease more 
strongly than the corresponding curve for the less compact EoS B series. As an example, while increasing the 
differential rotation parameter from A" 1 = 0.0 to A^ 1 = 0.43, the value of j3 c is lowered only by 6% for sequence B 
but by 30% for sequence II and even 39% for sequence A. Therefore, the onset of the CFS instability is also eased by 
the compactness of the neutron star when differential rotation is considered. 



2, Large Degree of Differential Rotation 



As we discussed in the previous section, equilibrium sequences with a high degree of differential rotation (A -1 > 1) 
terminate due to dynamical instabilities with respect to radial oscillations; see |49j for a study on the maximum mass 
of differentially rotating neutron stars. Therefore, it does not make sense to define the mass-shedding value /3 S , and 
correspondingly the same holds for the normalized quantity f3 c . Nevertheless, the sequences still reach the neutral 
point f3 c for the onset of the CFS instability of the fundamental quadrupolar mode. 

In Figure |8j we show f3 c of the quadrupolar f-mode for the three equations of state EoS B, EoS II and EoS A as 
well as for an additional one which we label EoS C and which has an intermediate polytropic exponent of T = 2.18 
and K = 400 in the same units as for the other equations of state. 

In contrast to the previous section, we continue to increase the degree of differential rotation up to A^ 1 = 1.75. 
When moving from rigid rotation to a small degree of differential rotation, the value of (3 C increases for each equation 
of state. Any sequence considered in this study reaches a local maximum which is located at A^ 1 = 0.89 for the 
B-sequence and A^ 1 = 0.78 for EoS C while EoS II and EoS A reach it considerably earlier at A^ 1 = 0.74 and 
A^ 1 — 0.69, respectively. 
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FIG. 8: The critical value /3 C of the 2 f-mode for four different equations of state. 

For higher degrees of differential rotation, the value of T/\W\ for the neutral point decreases. Again, similar to the 
observations we made for small degrees of differential rotation, this effect is more pronounced for the more compact 
models of the stiffer equations of state EoS II and EoS A. As one also can see from Figure [8j the critical value for 
the C, A and II sequences even drop clearly below the value for rigid rotation. For Eos A and EoS II, this happens 
roughly at around Ar 1 ss 1.2 while Ar 1 « 1.5 for EoS C. In contrast to this, while also decreasing in amplitude, the 
value of p c for the B-sequence does not drop below the corresponding value of the uniform rotation limit, at least not 
up to A^ 1 = 1.75 which is the endpoint of the sequences in Figure |8| 

A similar investigation has already been carried out for Newtonian polytropes [17] , and we qualitatively agree with 
their results. Further comparisons are not possible for various reasons; most crucial is the different definition of the 
parameter A which controls the degree of differential rotation. 

D. Rotational Modes in Differentially Rotating Stars 

Another very important class of oscillation modes are the r-modes, because they are generically CFS unstable 50J. 
Here, we only consider the I = m = 2 r-mode. Since the r-mode is purely axial in the non-rotating limit, it belongs 
to that class of modes whose pressure variation has odd parity under reflection with respect to the equatorial plane. 

In Figure [9] we show the frequencies of the r-mode along sequences using EoS B of for different degrees of differential 
rotation. We normalize them by the central angular velocity f2 c of the corresponding star, as it has already been done 
in [51] . This kind of normalization is suggested by the relation a — 4/3 which is valid in Newtonian theory, where 
a is the r-mode's angular frequency measured in the inertial frame and il the angular velocity of the star. 

First, we inspect the solid curve representing the rigidly rotation case, i.e. A^ 1 = 0. Concerning the value of 
a/Q c in the non-rotating limit, which is approximately 1.43, we observe a clear deviation from Newtonian theory. 
However, this is consistent with general relativistic corrections to the Newtonian result. It can be shown [52j that in 
the non-rotating limit the Newtonian value of 4/3 is increased by a term which is proportional to the frame-dragging 
potential uj. 

The curves corresponding to the differentially rotating sequences have lower values, because the central angular 
velocity increases compared to the uniformly rotating ones, and this effect becomes even stronger the higher the 
degree of differential rotation is. On the other hand, the angular velocity at the equator, O e , is lowered so that a 
normalization by this quantity would move the lines in the other direction. Finally, we chose to normalize by the 
central angular velocity in order to compare qualitatively with Newtonian findings [51] . We find that the qualitative 
behaviour is the same. 

Moreover, we explicitly computed the eigenfunction of the Su velocity perturbation. According to the Newtonian 
results in [51 , we find that the amplitude of the eigenfunction gets confined towards the surface of the star as 
differential rotation is increased. 

Earlier studies based on the slow rotation approximation also suggested that there might be a continuous spectrum 
for the r-modes [53H551 . More recently the same claim has been made in the study of differentially rotating relativistic 
stars in the slow-rotation approximation |19| . Here, in agreement with earlier results for uniformly rotating stars 
[38], 139] . we have found no trace of the continuum spectrum. This result does not exclude its existence, because 
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FIG. 9: Frequencies of the r-mode for several B sequences with different degrees of differential rotation. The frequencies are 
normalized by the central angular velocity fl c . 

the continuum can co-exist with discrete modes. Especially here we excite the oscillation by using an approximate 
eigenfunction of the mode under study, then by repeated mode recycling we force the star to oscillate on this specific 
mode alone in order to extract both the eigenfrequency and the eigenfunction with highest possible accuracy. This 
way we exclude excitation of any other mode or continuum, still we should stress that no numerical calculation of 
fast rotating relativistic star modes to date both in linear and non-linear form has shown the presence of a continuous 
spectrum. 



In this work we performed time evolutions of oscillations of differentially and rapidly rotating, relativistic stars using 
the Cowling approximation, and we were able to investigate both axisymmetric and non-axisymmetric perturbations. 
We extracted eigenfrequencies and eigenfunctions by applying Fourier transforms on the obtained time-series. The 
agreement with literature values is excellent for both non-axisymmetric perturbations of uniformly rotating stars and 
axisymmetric perturbations of differentially rotating stars. 

Our simulations show that a small degree of differential rotation causes the ratio of the critical T/|W|, where the 
f-mode becomes CFS unstable, to the value of T/|W| at the mass-sheeding to decrease. This effect is visible for each 
considered equation of state and is even stronger for the stiffer ones. It implies that neutron stars can be destabilized 
by the CFS instability even if their spinning frequency is rather low compared to the Kepler frequency. Furthermore, 
when moving on to a high degree of differential rotation, the critical curve of any considered sequence exhibits a 
local maximum which is in agreement with Newtonian findings. Afterwards, the critical value decreases and for the 
sequences based on the stiffer equations of state it decreases even below the corresponding value for rigid rotation; i. e. 
the CFS instability may act even earlier than in rigidly rotating stars as long as the star rotates with a sufficiently 
high degree of differential rotation. 

Summarized, we showed that both a small and a high degree of differential rotation plays an important role in the 
evolution of a neutron star. As explained in the introduction, neutron stars rotate differentially most probably directly 
after their birth, but may do so also at later times. Whenever differential rotation is present, the CFS instability 
is a promising candidate for the destabilization of a neutron star which may finally lead to a detectable signal of 
gravitational waves. 

Finally, we performed simulations on r-modes. Our results show a good qualitative agreement to Newtonian results. 
As in earlier studies which do not rely on the slow rotation approximation, we do not see any trace of a continuous 
spectrum. 

A very important issue related to differential rotation is the freedom to construct neutron star models with much 
higher values of T/| W|. This means that nearly all available equations of state can admit models which can be CFS- 
unstable which is not the case for uniformly rotating stars. This enhances the probability of having the CFS-instability 
working in the early phases of the proto-neutron star creation which will be a significant boost in our efforts to detect 
gravitational waves from isolated neutron stars and via such observations to infer and constrain fundamental neutron 
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star parameters. 
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